/*****************************   model.h   **********************************
* Author:        Gena Tang
* Date created:  2012
* Description:
* This header file contains class declarations and other definitions for the 
* model class library of GeneMark model
******************************************************************************/

#ifndef MODEL_H
#define MODEL_H


#include <string>
#include <vector>
#include <map>
using namespace std;



enum NUCLEOTIDE { A, C, G, T};



class ModelFile {                // GeneMark .mod file reader

public:
        ModelFile();         // Constructor
        void read(string modfile); //default order is 2
		void AbsToFirst( vector<double> & abs, vector<double>& relative );
		void AbsToLast( vector<double> & abs);
		void getCodLengthDistr();
		void getNonLengthDistr();
        //void read(string modfile, int order);         

        int order; //model order
		int msize; //parameter size m= 4, 16, 64...
        int dur_non; // noncoding duration
		double codDecay;
		double noncodDecay;
		double codToCod;
        vector<double> P_1; //native coding model frame 1 
		vector<double> P_2; //native coding model frame 2 
		vector<double> P_3; //native coding model frame 3 
        vector<double> H_1; //heuristic coding model frame 1
		vector<double> H_2; //heuristic coding model frame 2
		vector<double> H_3; //heuristic coding model frame 3
        vector<double> Q; //noncoding model
		vector<double> P; //coding model no frames. (avg of P_1, P_2, P_3).

		//relative probablity. P[aaa]+P[aac]+P[aag]+P[aat] = 1 not used!!!
		vector<double> P_1_rel; //native coding model frame 1 
		vector<double> P_2_rel; //native coding model frame 2 
		vector<double> P_3_rel; //native coding model frame 3 
        vector<double> H_1_rel; //heuristic coding model frame 1
		vector<double> H_2_rel; //heuristic coding model frame 2
		vector<double> H_3_rel; //heuristic coding model frame 3
        vector<double> Q_rel; //noncoding model
		vector<double> codLengthDistr;
		vector<double> nonLengthDistr;

        vector<double> P_start; //start codon usage
        vector<double> Q_4; // zero order model for non-coding. Q_4[0] is p(A);

		vector<string> stop_codons;
		vector<string> start_codons;

		int minc; //min length of coding
		int maxc;
		int maxn;
        int spacer_l; //spacer length 
        int motif_l; //RBS site length
        vector<double> spacer; //spacer distribution, spacer[0] is closest to gene start
        vector< vector <double> > RBS; //motif matrix, RBS[0][0] is furthest to gene start
		vector< vector <double> > StartM; //start codon

private:
        double durationFromModel( vector<double>& arr );
        void removeStopCodons( vector<double>& arr );
        void normArr( vector<double>& arr );
        void getZero( vector<double>& arr, vector<double>& outarr ); //get zero order model from 2 order.
}; 


class Sequence {                // .fna/.fasta file reader

public:
        Sequence();         // null Constructor
        void read(string fastafile); //read sequence, only one seuqence can be read
        string name;
        string seq; 

private:

}; 


class Gene {                

public:
        Gene(); // Constructor
        Gene(string line);         // Constructor, input is one line of .lst annotation

        void set_gene_seq(string & seq);
        void set_upstream_seq(string & seq, int l);

        double get_codon_score(double * P_Q);
        double get_start_score(double * P_Q);
        double get_length_score(vector<double>& P_Q_length_logodd);


        int id;
        char strand;
        int gene_left;
        int gene_right;
        int gene_length;
        int gene_class;
        string gene_seq; // seq represented by 0-3
        string upstream_seq; // seq represented by 0-3
        double codon_score;
        double length_score;
        double motif_score;
        double score;
        double p_value;

private:
        void seqToInt(string & seq);
}; 





void printVector(string name, vector<double>& arr);
string reverse(string s);
string revcomp(string s);
string revcompToNT(string str);
double* VectorToArray( vector<double>& arr );
int* VectorToArray( vector<int>& arr );
double* LogOdd( double* F, double* S, int size );
void LengthLogOdd( vector<double>& logodd, int min, int max, double gamma, double exp );
double NormExpConstant( double betta, int min, int max );
double NormGammConstant( double betta, int min, int max );
double GetExp( int i, double betta );
double GetGamma( int i, double betta );
void avg(vector<double>& p, vector<double>& p1, vector<double>& p2, vector<double>& p3);

#endif